// REPLICATON FILE FOR "ASSET RETURNS AND LABOR SUPPLY IN A PRODUCTION ECONOMY", IVAN JACCARD, APRIL 2013
// ARTICLE ACCEPTED FOR PUBLICATION IN THE JOURNAL OF MONEY, CREDIT AND BANKING
// REQUIRES DYNARE 4.3.3

//---------------------------------------------------------------------
// 1. Variable declaration
//---------------------------------------------------------------------

var c, n, k, x, lambda, phi, y, i, A, rf, q, d, pe, w, comp, epann, rfann, realret_e, gy, gc, gi, gn, gw, logcomp, logc, logn, ep;
varexo ez ;                          


//---------------------------------------------------------------------
// 2. Parameter declaration and calibration
//---------------------------------------------------------------------

parameters wlambda, WELAST, beta2, alpha, delta, gamma, sigma, rho, beta, v, m, th1, th2, eps, kyratio, Si, bstar, psi, N, L, Sc, y_sts, c_sts, i_sts, k_sts, x_sts, phi_sts, lambda_sts, d_sts, pe_sts ;

//**********************
rho = 0.999      ;
bstar = 0.98477  ; 
eps =3.8526      ;  
m= 0.50          ; 
//***********************

alpha      = 0.36;        // capital share in production function         
delta      = 0.025;      // quarterly depreciation rate       
gamma      = 1.007;       //quarterly growth rate 
sigma      = 1.0;         // curvature parameter    

th1        = (gamma-1+delta)^(eps);            // adjustment cost param1          
th2        = -(gamma-1+delta)*(eps/(1-eps)) ;  // adjustment cost param2
beta       = bstar*gamma^(-sigma);             // modified discount factor 1, assets that are growing in the steady state
beta2      = bstar*gamma^(1-sigma);            // modified discount factor 2, assets that are not growing in the steady state
v          = 3.7;                              // parameter governing Frish elasticity, set to have Frisch close but below 1
N          = 0.2;                              // steady state hours worked
L          = 1-0.2;                            // steady state leisure


// STEADY STATE IMPLICATONS

Sc         = 1-(beta*alpha*(gamma-(1-delta)))/(1-beta*(1-delta));      // consumption to output ratio
psi        = (N/(1-alpha))*Sc*v*(L^(v-1))-(L^v);                       // labor supply parameter
kyratio    = beta*alpha/(1-beta*(1-delta));
Si         = 1-Sc;

// CHECK STEADY STATE
y_sts      = ((kyratio)^(alpha/(1-alpha)))*N;
c_sts      = Sc*y_sts;
i_sts      = Si*y_sts;
k_sts      = kyratio*y_sts;
x_sts      = ((1-m)/(gamma-m))*(c_sts*(psi+L^v));
phi_sts    = -beta*(((c_sts*(psi+L^v)-x_sts)^-sigma)/(1-beta*m));
lambda_sts = ((c_sts*(psi+((1-N)^v))-x_sts)^-sigma)*(psi+((1-N)^v))+(1-m)*phi_sts*(((1-N)^v)+psi);
d_sts      =  alpha*y_sts-i_sts;
pe_sts     =  (beta2/(1-beta2))*d_sts;   

//FRISCH
RHO    = 1-(1-m)/(gamma-m);
GAMMA  = 1-(1-m)*beta/(1-beta*m);
eC     = 1-sigma/(RHO*GAMMA);
eN     = ((v-1)-v*((L^v)/(psi+(L^v)))*(sigma/(RHO*GAMMA)))*(N/(1-N));
wN     = v*((L^v)/(psi+(L^v)))*(1-(RHO*GAMMA)/sigma)*(N/(1-N));
wlambda = RHO*GAMMA/sigma;
FRISCH = 1/(eC*wN-eN);
WELAST = (1+eC*wlambda)/(eC*wN-eN);

//---------------------------------------------------------------------
// 3. Model declaration
//---------------------------------------------------------------------

model;  

((c*(psi+((1-n)^v))-x(-1))^-sigma)*(psi+((1-n)^v))+(1-m)*phi*(((1-n)^v)+psi)     = lambda;

((c*(((1-n)^v)+psi)-x(-1))^-sigma)*v*c*((1-n)^(v-1))+(1-m)*phi*v*c*((1-n)^(v-1)) = lambda*(1-alpha)*y/n;

1            =   q*th1*(i/k(-1))^(-eps);

phi          =   m*beta*phi(+1)-beta*((c(+1)*(psi+((1-n(+1))^v))-x)^(-sigma));

lambda*q     =   beta*q(+1)*(lambda(+1))*((1-delta)+(th1/(1-eps))*((i(+1)/k)^(1-eps))+th2-th1*(i(+1)/k)^(1-eps))+beta*((lambda(+1))*(alpha*(y(+1)/k)));

gamma*k      =   (1-delta)*k(-1) + (k(-1))*((th1/(1-eps))*((i/k(-1))^(1-eps))+th2);

gamma*x      =   m*x(-1)+(1-m)*c*(psi+((1-n)^v));

y            =   c+i;

y            =   A*(n^(1-alpha))*(k(-1)^alpha);

log(A)       =   rho*log(A(-1))+ez;

w            =   (1-alpha)*A*(k(-1)^alpha)*(n^(-alpha));

comp         =   c*(psi+((1-n)^v));

// ASSET PRICING IMPLICATIONS

pe*lambda    =   beta2*lambda(+1)*(d(+1)+pe(+1));

1/rf         =   beta2*lambda(+1)/lambda;

d            =   alpha*y-i;

rfann        =  (rf-1)*400;

epann        = ((pe(+1)+d(+1))/pe-rf)*400;

realret_e    =  (pe+d)/pe(-1);

ep           =   ((pe(+1)+d(+1))/pe-rf);

// GROWTH RATES

gy          =  (log(y)-log(y(-1)))*100;

gc          =  (log(c)-log(c(-1)))*100;

gi          =  (log(i)-log(i(-1)))*100;

gn          =  (log(n)-log(n(-1)))*100;

gw          =  (log(w)-log(w(-1)))*100;

logcomp     =  log(comp);

logc        =  log(c); 

logn        = log(n);

end;

//---------------------------------------------------------------------
// 4. Initial values and steady state
//---------------------------------------------------------------------

initval;                      

n       = N;
A       = 1;
y       = ((kyratio)^(alpha/(1-alpha)))*n;
c       = Sc*y;
k       = kyratio*y;
x       = x_sts;
phi     = phi_sts;
lambda  = lambda_sts;
rf      = 1/beta2;
q       = 1;
i       = i_sts;
d       = d_sts;
pe      = pe_sts;
w       = (1-alpha)*A*(k^alpha)*(n^(-alpha));
comp    = c*(psi+((1-n)^v));
rfann   = (rf-1)*400;
realret_e = (pe+d)/pe;
epann   = 0;
gy      = 0;
gc      = 0;
gi      = 0;
gn      = 0;
gw      = 0;
logcomp =  log(comp);
logc    =  log(c); 
logn    = log(n);
ep      = 0;
end;

//---------------------------------------------------------------------
// 5. Shock declaration  
//                       
//---------------------------------------------------------------------

shocks;
var ez; stderr 0.0072;  
end;

stoch_simul (order=2, irf=0) epann, gy, gc, gi, gn, gw;


